knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, fig.height = 4, fig.width = 7)
library(terra)
library(sf)
library(oharac)
library(data.table)
library(tidyverse)
library(biscale) ### for bivariate legend
library(cowplot)
library(here)
source(here('common_fxns.R'))Create Fig. 2 for manuscript: A) a panel of cumulative impacts, B) & C) separated into climate and non-climate, D) quartile overlaps
Results from scripts in folder
2_process_spp_vuln_and_impacts
Read in the output rasters for species-level impacts (unweighted) - total, climate, and non-climate; reclassify based on quartile/quintile for biplot.
ocean_map <- rast(here('_spatial/ocean_area_mol.tif'))
ocean_df <- as.data.frame(ocean_map, xy = TRUE)
land_sf <- read_sf(here('_spatial/ne_10m_land/ne_10m_land_no_casp.shp')) %>%
st_transform(st_crs(ocean_map))
nspp_mask <- rast(here('_output/nspp_maps/species_richness.tif'))
nspp_vals <- values(nspp_mask)
values(nspp_mask)[nspp_vals < 20] <- NA
chi_gps_spp <- list.files(here('_output/cumulative_impact_maps'),
pattern = 'species_group_chi',
full.names = TRUE)
climate_map_spp <- rast(chi_gps_spp[1])
nonclimate_map_spp <- rast(chi_gps_spp[2:4]) %>%
sum(na.rm = TRUE)
chi_total_map_spp <- rast(here('_output/cumulative_impact_maps/chi_species.tif'))Squashing extreme outliers!
squash_outliers <- quantile(values(chi_total_map_spp), .9999, na.rm = TRUE)
cc_spp_df <- as.data.frame(climate_map_spp, xy = TRUE) %>%
setNames(c('x', 'y', 'cc_spp')) %>%
mutate(cc_spp = ifelse(cc_spp < squash_outliers, cc_spp, squash_outliers))
noncc_spp_df <- as.data.frame(nonclimate_map_spp, xy = TRUE) %>%
setNames(c('x', 'y', 'noncc_spp')) %>%
mutate(cc_spp = ifelse(noncc_spp < squash_outliers, noncc_spp, squash_outliers))
chi_spp_df <- as.data.frame(chi_total_map_spp, xy = TRUE) %>%
setNames(c('x', 'y', 'chi')) %>%
mutate(chi = ifelse(chi < squash_outliers, chi, squash_outliers))
chi_lims <- c(0, round(max(chi_spp_df$chi), 2))
chi_brks <- c(0, .25, .5, 0.75, round(max(chi_spp_df$chi), 2))
chi_lbls <- sprintf('%.2f', chi_brks)
chi_lbls[length(chi_lbls)] <- paste0('≥ ', chi_lbls[length(chi_lbls)])
xfm <- function(x) {x^1} ### function to warp CHI scale for visual impactglobe_bbox <- rbind(c(-180, -90), c(-180, 90),
c(180, 90), c(180, -90), c(-180, -90))
globe_border <- st_polygon(list(globe_bbox)) %>%
st_sfc(crs = 4326) %>%
st_sf(data.frame(rgn = 'globe', geom = .)) %>%
smoothr::densify(max_distance = 0.5) %>%
st_transform(crs = crs(ocean_map))inset_bbox_df <- tribble(
~xmin, ~xmax, ~ymin, ~ymax, ~loc,
-1.0e6, +3.5e6, +3.5e6, +8.0e6, 'europe',
+8.5e6, +13.5e6, -0.0e6, +5.0e6, 'se_asia',
-9.5e6, -5.5e6, +0.5e6, +4.5e6, 'caribbean'
)
panel_a_chi <- ggplot() +
### plot data:
geom_raster(data = chi_spp_df, aes(x, y, fill = xfm(chi))) +
### continents:
geom_sf(data = land_sf,
fill = 'grey96', color = NA,
size = .1) +
geom_sf(data = globe_border,
fill = NA, color = 'grey70',
size = .1) +
geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
color = 'red', fill = NA, linewidth = 0.2) +
scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_lbls) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme_void() +
theme(axis.text = element_blank(),
axis.title = element_blank(),
legend.text = element_text(hjust = 0.5, vjust = 1, size = 7, angle = 0),
legend.title = element_blank(),
legend.key.height = unit(.16, 'in'),
legend.key.width = unit(.32, 'in'),
legend.direction = 'horizontal')
# panel_a_chi
# asdf <- get_legend(panel_a_chi)plot_inset <- function(chi_df, bbox_df) {
### for chi_df, rename the chi column to just 'chi'
### for bbox_df, filter to the appropriate row
ggplot() +
### plot data:
geom_raster(data = chi_df, aes(x, y, fill = xfm(chi)), show.legend = FALSE) +
scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_lbls) +
scale_x_continuous(expand = c(0, 0), limits = c(bbox_df$xmin[1], bbox_df$xmax[1])) +
scale_y_continuous(expand = c(0, 0), limits = c(bbox_df$ymin[1], bbox_df$ymax[1])) +
theme_void() +
theme(panel.border = element_rect(color = 'black', fill = NA, linewidth = .5),
panel.background = element_rect(color = NA, fill = 'grey96'))
}
inset_a_i <- plot_inset(chi_df = chi_spp_df, bbox_df = inset_bbox_df %>% filter(loc == 'europe'))
inset_a_ii <- plot_inset(chi_df = chi_spp_df, bbox_df = inset_bbox_df %>% filter(loc == 'caribbean'))
inset_a_iii <- plot_inset(chi_df = chi_spp_df, bbox_df = inset_bbox_df %>% filter(loc == 'se_asia'))panel_b_cc <- ggplot() +
### plot data:
geom_raster(data = cc_spp_df, aes(x, y, fill = xfm(cc_spp))) +
### continents:
geom_sf(data = land_sf,
fill = 'grey96', color = NA,
size = .1) +
geom_sf(data = globe_border,
fill = NA, color = 'grey70',
size = .1) +
# geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
# color = 'red', fill = NA, linewidth = 0.2) +
scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_brks) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme_void() +
theme(axis.text = element_blank(),
axis.title = element_blank())
# panel_b_ccpanel_c_noncc <- ggplot() +
### plot data:
geom_raster(data = noncc_spp_df, aes(x, y, fill = xfm(noncc_spp))) +
### continents:
geom_sf(data = land_sf,
fill = 'grey96', color = NA,
size = .1) +
geom_sf(data = globe_border,
fill = NA, color = 'grey70',
size = .1) +
# geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
# color = 'red', fill = NA, linewidth = 0.2) +
scale_fill_viridis_c(limits = xfm(chi_lims), breaks = xfm(chi_brks), labels = chi_brks) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme_void() +
theme(axis.text = element_blank(),
axis.title = element_blank())
# panel_c_nonccUse a four-level bivariate color scale to capture spatial overlap of quartiles for climate and non-climate stressors.
Reclassify values using ntile to get quartile
values:
cc_qtile_spp <- noncc_qtile_spp <- climate_map_spp ### initialize two new maps
values(cc_qtile_spp) <- ntile(values(climate_map_spp), 4)
values(noncc_qtile_spp) <- ntile(values(nonclimate_map_spp), 4)cc_spp_qtile_df <- as.data.frame(cc_qtile_spp, xy = TRUE) %>%
setNames(c('x', 'y', 'cc_spp'))
noncc_spp_qtile_df <- as.data.frame(noncc_qtile_spp, xy = TRUE) %>%
setNames(c('x', 'y', 'noncc_spp'))
cc_noncc_spp_qtile_df <- full_join(cc_spp_qtile_df, noncc_spp_qtile_df, by = c('x', 'y')) %>%
mutate(bi_class = paste(cc_spp, noncc_spp, sep = '-'))
cc_noncc_summary <- cc_noncc_spp_qtile_df %>%
group_by(cc_spp, noncc_spp) %>%
summarize(n = n(),
pct = n() / nrow(cc_noncc_spp_qtile_df),
.groups = 'drop') %>%
mutate(pct = round(pct * 100, 1),
lbl = paste0(pct, '%'))
knitr::kable(cc_noncc_summary, digits = 3)| cc_spp | noncc_spp | n | pct | lbl |
|---|---|---|---|---|
| 1 | 1 | 324581 | 8.8 | 8.8% |
| 1 | 2 | 300192 | 8.2 | 8.2% |
| 1 | 3 | 179527 | 4.9 | 4.9% |
| 1 | 4 | 115698 | 3.1 | 3.1% |
| 2 | 1 | 279191 | 7.6 | 7.6% |
| 2 | 2 | 250271 | 6.8 | 6.8% |
| 2 | 3 | 231144 | 6.3 | 6.3% |
| 2 | 4 | 159392 | 4.3 | 4.3% |
| 3 | 1 | 211177 | 5.7 | 5.7% |
| 3 | 2 | 178541 | 4.9 | 4.9% |
| 3 | 3 | 264828 | 7.2 | 7.2% |
| 3 | 4 | 265452 | 7.2 | 7.2% |
| 4 | 1 | 105049 | 2.9 | 2.9% |
| 4 | 2 | 190994 | 5.2 | 5.2% |
| 4 | 3 | 244499 | 6.6 | 6.6% |
| 4 | 4 | 379456 | 10.3 | 10.3% |
For the bivariate color plot, one axis of colors represents climate, the other represents non-climate; for corner values (max agreement and max disagreement), use diverging color palettes; for intermediate values use a light shading to indicate variation more subtly.
map_pal_raw <- bi_pal(pal = 'PurpleOr', dim = 4, preview = FALSE)
map_pal_mtx <- matrix(map_pal_raw, nrow = 4, ncol = 4)
map_pal_mtx[3, ] <- colorspace::lighten(map_pal_mtx[3, ], .1)
map_pal_mtx[2, ] <- colorspace::lighten(map_pal_mtx[2, ], .2)
map_pal_mtx[1, ] <- colorspace::lighten(map_pal_mtx[1, ], .3)
map_pal_mtx[ , 3] <- colorspace::lighten(map_pal_mtx[ , 3], .1)
map_pal_mtx[ , 2] <- colorspace::lighten(map_pal_mtx[ , 2], .2)
map_pal_mtx[ , 1] <- colorspace::lighten(map_pal_mtx[ , 1], .3)
map_pal_mtx[1, 1] <- '#ffffee'
map_pal <- as.vector(map_pal_mtx) %>% setNames(names(map_pal_raw))
panel_lgd_raw <- bi_legend(pal = map_pal, dim = 4,
x = 'Climate',
y = 'Non-climate')
si_fig <- panel_lgd_raw +
geom_text(data = cc_noncc_summary %>%
filter(cc_spp + noncc_spp < 6),
aes(x = cc_spp, y = noncc_spp, label = lbl),
hjust = .5, vjust = .5, color = 'black', size = 3) +
geom_text(data = cc_noncc_summary %>%
filter(cc_spp + noncc_spp >= 6),
aes(x = cc_spp, y = noncc_spp, label = lbl),
hjust = .5, vjust = .5, color = 'white', size = 3) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme(axis.title = element_text(hjust = 0, size = 12, face = 'bold', family = 'Arial'),
panel.border = element_rect(color = 'black', fill = NA))
si_figggsave(filename = 'figS3A_spp_legend_w_pcts.png', plot = si_fig, height = 3, width = 3, units = 'in', dpi = 300)panel_d_biplot <- ggplot() +
### plot data:
geom_raster(data = cc_noncc_spp_qtile_df, aes(x, y, fill = bi_class), show.legend = FALSE) +
### continents:
geom_sf(data = land_sf,
fill = 'grey96', color = NA,
size = .1) +
geom_sf(data = globe_border,
fill = NA, color = 'grey70',
size = .1) +
# geom_rect(data = inset_bbox_df, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
# color = 'red', fill = NA, linewidth = 0.2) +
bi_scale_fill(pal = map_pal, dim = 4) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme_void() +
theme(axis.text = element_blank(),
axis.title = element_blank())
# panel_d_biplotfig2_file <- here('5_ms_figs/fig2_mapping_impacts_spp.png')
panel_a_map <- get_panel(panel_a_chi)
panel_b_map <- get_panel(panel_b_cc)
panel_c_map <- get_panel(panel_c_noncc)
panel_d_map <- get_panel(panel_d_biplot)
panel_d_lgd <- get_panel(panel_lgd_raw)
panel_abc_lgd <- get_legend(panel_a_chi)
panel_e_inset <- get_panel(inset_a_i)
panel_f_inset <- get_panel(inset_a_ii)
panel_g_inset <- get_panel(inset_a_iii) ### parameterize panels for easier changing
fig_h <- 8.5
fig_w <- 6.25
ar <- fig_h / fig_w
p_h <- 0.24
p_w <- p_h * 2 * ar
i_h <- p_h * 0.8
i_w <- i_h * ar
a <- c(x = 0.00, y = 0.75)
b <- c(x = 0.00, y = 0.50)
c <- c(x = 0.00, y = 0.25)
d <- c(x = 0.00, y = 0.0)
e <- c(x = .968 - i_w, y = 0.65)
f <- c(x = .968 - i_w, y = 0.45)
g <- c(x = .968 - i_w, y = 0.25)
fig2_panels <- ggdraw() +
draw_plot(panel_a_map, x = a['x'], y = a['y'], height = p_h, width = p_w) +
draw_plot(panel_b_map, x = b['x'], y = b['y'], height = p_h, width = p_w) +
draw_plot(panel_c_map, x = c['x'], y = c['y'], height = p_h, width = p_w) +
draw_plot(panel_d_map, x = d['x'], y = d['y'], height = p_h, width = p_w) +
draw_plot(panel_e_inset, x = e['x'], y = e['y'], height = i_h, width = i_w) +
draw_plot(panel_f_inset, x = f['x'], y = f['y'], height = i_h, width = i_w) +
draw_plot(panel_g_inset, x = g['x'], y = g['y'], height = i_h, width = i_w)
fig2_labels <- fig2_panels +
draw_label('A.', x = a['x'], y = a['y'] + p_h - .05, hjust = 0, vjust = 1,
size = 12, fontfamily = 'Arial', fontface = 'bold') +
draw_label('B.', x = b['x'], y = b['y'] + p_h - .05, hjust = 0, vjust = 1,
size = 12, fontfamily = 'Arial', fontface = 'bold') +
draw_label('C.', x = c['x'], y = c['y'] + p_h - .05, hjust = 0, vjust = 1,
size = 12, fontfamily = 'Arial', fontface = 'bold') +
draw_label('D.', x = d['x'], y = d['y'] + p_h - .05, hjust = 0, vjust = 1,
size = 12, fontfamily = 'Arial', fontface = 'bold') +
draw_label('E.', x = e['x'] - .01, y = e['y'] + i_h, hjust = 1, vjust = 1,
size = 12, fontfamily = 'Arial', fontface = 'bold') +
draw_label('F.', x = f['x'] - .01, y = f['y'] + i_h, hjust = 1, vjust = 1,
size = 12, fontfamily = 'Arial', fontface = 'bold') +
draw_label('G.', x = g['x'] - .01, y = g['y'] + i_h, hjust = 1, vjust = 1,
size = 12, fontfamily = 'Arial', fontface = 'bold')### parameterize legends for easier changing of label positions:
l1 <- c(x = p_w + .02, y = 0.85, h = .10, w = .97 - p_w)
l2 <- c(x = p_w + .048, y = 0.06, h = .12, w = .12 * ar)
fig2_legends <- fig2_labels +
### draw legends:
draw_plot(panel_abc_lgd, x = l1['x'], y = l1['y'],
height = l1['h'], width = l1['w']) +
draw_plot(panel_d_lgd, x = l2['x'], y = l2['y'],
height = l2['h'], width = l2['w']) +
### Legend label for CHI:
draw_label('Cumulative human impacts\n(A-C, E-G)',
x = l1['x'] + l1['w']/2, y = l1['y'] + l1['h'],
hjust = 0.5, vjust = 0.5, angle = 0,
size = 10, color = 'black', fontfamily = 'Arial', fontface = 'bold') +
### Legend labels for biscale:
draw_label('Impact\nquartiles (D)',
x = l2['x'] + .01, y = l2['y'] + l2['h'],
hjust = 0, vjust = 0, angle = 0,
size = 10, color = 'black', fontfamily = 'Arial', fontface = 'bold') +
draw_label('Climate →',
x = l2['x'], y = l2['y'] - .005,
hjust = 0, vjust = 1, angle = 0,
size = 8, color = 'black', fontfamily = 'Arial', fontface = 'bold') +
draw_label('Non-climate →',
x = l2['x'] - .005, y = l2['y'],
hjust = 0, vjust = 0, angle = 90,
size = 8, color = 'black', fontfamily = 'Arial', fontface = 'bold')
ggsave(fig2_file, width = fig_w, height = fig_h, units = 'in', dpi = 300)
ggsave(fig2_file %>% str_replace('.png$', '.tif'), width = fig_w, height = fig_h, units = 'in', dpi = 300)
knitr::include_graphics(fig2_file)